The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file
Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here
#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.3 ✔ forcats 0.5.2 ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))
processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd
|
| | 0%
|
|............. | 5%
|
|........................... | 10%
|
|........................................ | 14%
|
|..................................................... | 19%
|
|.................................................................. | 24%
|
|................................................................................ | 29%
|
|............................................................................................. | 33%
|
|.......................................................................................................... | 38%
|
|........................................................................................................................ | 43%
|
|..................................................................................................................................... | 48%
|
|.................................................................................................................................................. | 52%
|
|............................................................................................................................................................... | 57%
|
|............................................................................................................................................................................. | 62%
|
|.......................................................................................................................................................................................... | 67%
|
|....................................................................................................................................................................................................... | 71%
|
|..................................................................................................................................................................................................................... | 76%
|
|.................................................................................................................................................................................................................................. | 81%
|
|............................................................................................................................................................................................................................................... | 86%
|
|............................................................................................................................................................................................................................................................ | 90%
|
|.......................................................................................................................................................................................................................................................................... | 95%
|
|.......................................................................................................................................................................................................................................................................................| 100%
output file: /tmp/RtmpPXtcTe/file315a4f0aaec
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
Attaching package: ‘flextable’
The following object is masked from ‘package:purrr’:
compose
Attaching package: ‘ftExtra’
The following object is masked from ‘package:flextable’:
separate_header
Warning: Assuming taxa are rows
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3
library(cowplot)
library(flextable)
library(ftExtra)
This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.
Reshape back into an ASV matrix, but this time correcting for total abundance
nonSpikes %>% head
raDf <- nonSpikes %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")
raMat1 <- raMat %>% as.matrix()
countMat <- nonSpikes %>%
pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
column_to_rownames("ID") %>% as.matrix()
seqDep <- countMat %>% apply(1, sum)
min(seqDep)
[1] 852
sampleRichness <- rarefy(countMat, min(seqDep))
rarefy everything to the minimum depth (852)
countRare <- rrarefy(countMat, min(seqDep))
Gamma diversity
specpool(countRare)
Doesn’t finish
#specpool(countMat)
All richness estimates
richnessRare <- estimateR(countRare)
Shannon diversity
shan <- diversity(countRare)
shan
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53
4.382762 5.084370 4.678155 5.926793 5.237899 3.867372 5.596986 4.563925 4.967078 4.765325 5.298189 4.815413 4.407171 4.665767 4.535506 5.052774 5.446895 5.077076 4.960559 3.757001 4.947850 4.690475 4.971215 5.197389 4.867073 4.144527
3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20
4.317769 4.955861 3.441315 5.775953 5.283504 5.340186 5.501008 5.042675 5.026191 4.920511 4.302272 4.376545 4.946276 4.475293 4.367380 4.652247 4.290391 4.321040 5.034992 4.637150 5.244157 4.021406 4.520193 2.805419 4.545389 4.648545
4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
4.747469 4.443062 4.354328 4.570223 4.145389 4.095409 4.085605 4.167853 4.681186 5.103463 5.441884 5.083902 4.899043 4.294475 4.887705 4.875368 4.208047 4.295552 4.779480 4.864873 5.423950 4.773348 5.107284
Evenness
pielouJ <- shan/richnessRare["S.chao1",]
pielouJ
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5
0.011295779 0.008105812 0.012287837 0.003007367 0.005702428 0.042036654 0.008402445 0.010684687 0.007005046 0.008998109 0.006977728 0.009101954 0.014352266 0.009499526 0.014455797 0.006241515 0.007677090 0.008637849 0.012140483 0.017485268 0.008113446 0.010226522 0.006023986 0.009093441
3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500
0.009069482 0.014334723 0.012184340 0.009141746 0.064686366 0.004660697 0.007687693 0.005758500 0.006703955 0.008822806 0.008114157 0.009782327 0.012208135 0.010270972 0.008396074 0.009961698 0.009821496 0.008589214 0.008862728 0.011374560 0.008911489 0.009510012 0.007729699 0.011732768
4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2
0.009779128 0.124685269 0.012987829 0.009992095 0.013958218 0.014609567 0.010345631 0.012334306 0.014569903 0.017536336 0.012430848 0.015666936 0.007652523 0.008220342 0.005410565 0.009486091 0.009221795 0.016242008 0.011576142 0.009351730 0.015098839 0.009240293 0.007351504 0.005571680
C_5P1B_20 C_5P1B_500 C_5P1B_53
0.003605969 0.007392898 0.005904266
diversityData <- sampleData %>%
left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
arrange(Size_Class)
Parameters for all plots
plotSpecs <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
Observed species counts, on rarefied data
plotObs <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs
Estemated Chao1 Richness
plotChao1 <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) +
scale_y_log10() +
ylab("Richness (Chao1)")
plotChao1
Shannon diversity
plotShan <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
ylab("Diversity (Shannon H)") +
lims(y = c(2.5, 6))
plotShan
Evenness
plotPielou <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou
All plots together
plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha
ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)
Rarefied observed species numbers
obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)
Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-224.389 -40.954 2.457 53.730 199.561
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 292239.870 138034.341 2.117 0.037854 *
log(Size_Class) 32.107 8.363 3.839 0.000271 ***
I(log(Size_Class)^2) -6.364 1.556 -4.090 0.000115 ***
lat -15220.663 7192.791 -2.116 0.037947 *
I(lat^2) 198.241 93.636 2.117 0.037856 *
depth 5.152 3.312 1.555 0.124459
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 78.62 on 69 degrees of freedom
Multiple R-squared: 0.2554, Adjusted R-squared: 0.2015
F-statistic: 4.734 on 5 and 69 DF, p-value: 0.0008913
Rarified chao1 estimates
chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-539.32 -131.99 -38.59 128.05 1210.56
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.071e+06 4.831e+05 2.217 0.02989 *
log(Size_Class) 8.568e+01 2.927e+01 2.927 0.00463 **
I(log(Size_Class)^2) -1.801e+01 5.446e+00 -3.307 0.00150 **
lat -5.582e+04 2.517e+04 -2.218 0.02988 *
I(lat^2) 7.272e+02 3.277e+02 2.219 0.02978 *
depth 2.227e+01 1.159e+01 1.921 0.05884 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 275.2 on 69 degrees of freedom
Multiple R-squared: 0.1982, Adjusted R-squared: 0.1401
F-statistic: 3.411 on 5 and 69 DF, p-value: 0.008206
As above but without latitude and depth
chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2),
data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-466.78 -164.77 -27.36 117.52 1308.99
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 568.437 50.075 11.352 < 2e-16 ***
log(Size_Class) 86.353 29.738 2.904 0.00489 **
I(log(Size_Class)^2) -18.426 5.532 -3.331 0.00137 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 279.9 on 72 degrees of freedom
Multiple R-squared: 0.1341, Adjusted R-squared: 0.1101
F-statistic: 5.576 on 2 and 72 DF, p-value: 0.005603
shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(shanMod)
Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-1.35466 -0.22941 0.03062 0.30917 0.76763
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.683e+03 8.051e+02 2.090 0.04028 *
log(Size_Class) 1.990e-01 4.878e-02 4.079 0.00012 ***
I(log(Size_Class)^2) -3.661e-02 9.075e-03 -4.034 0.00014 ***
lat -8.744e+01 4.195e+01 -2.084 0.04083 *
I(lat^2) 1.139e+00 5.461e-01 2.085 0.04079 *
depth 2.099e-02 1.932e-02 1.087 0.28095
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4586 on 69 degrees of freedom
Multiple R-squared: 0.2863, Adjusted R-squared: 0.2346
F-statistic: 5.537 on 5 and 69 DF, p-value: 0.000241
pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)
Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-0.014482 -0.005102 -0.002446 0.000812 0.099932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.506e+01 2.657e+01 -0.943 0.3490
log(Size_Class) -3.808e-03 1.610e-03 -2.365 0.0208 *
I(log(Size_Class)^2) 6.620e-04 2.996e-04 2.210 0.0304 *
lat 1.305e+00 1.385e+00 0.942 0.3493
I(lat^2) -1.697e-02 1.803e-02 -0.941 0.3498
depth -3.651e-04 6.377e-04 -0.573 0.5688
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01514 on 69 degrees of freedom
Multiple R-squared: 0.09266, Adjusted R-squared: 0.02691
F-statistic: 1.409 on 5 and 69 DF, p-value: 0.2318
uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2
predict(obsMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
231.2129 231.2129 197.7419 197.7419 213.1276 161.0777 161.0777 193.7155 203.6695 305.0142 305.0142 271.5432 271.5432 286.9290 234.8790 234.8790 267.5169 267.5169 334.5602 334.5602 301.0891 301.0891 316.4749 264.4250 264.4250 297.0628 307.0168 307.0168 338.4386 338.4386 304.9676 304.9676
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
320.3533 320.3533 268.3034 268.3034 300.9413 300.9413 326.5221 293.0510 293.0510 308.4368 308.4368 256.3869 256.3869 256.3869 289.0247 289.0247 298.9787 298.9787 294.4742 294.4742 261.0032 261.0032 276.3890 276.3890 224.3390 224.3390 224.3390 256.9769 256.9769 266.9308 266.9308 253.1023
65 66 67 68 69 70 71 72 73 74 75
219.6313 219.6313 235.0171 235.0171 182.9671 182.9671 182.9671 215.6050 215.6050 225.5589 225.5589
$se.fit
[1] 27.50522 27.50522 26.75531 26.75531 26.61253 28.31263 28.31263 29.90239 34.68922 19.29454 19.29454 18.66638 18.66638 17.52990 20.42368 20.42368 21.91454 21.91454 19.62244 19.62244 19.15178 19.15178 17.56592 20.47128 20.47128 21.61090 28.71618 28.71618 19.62022 19.62022 19.11805
[32] 19.11805 17.28640 17.28640 19.98620 19.98620 20.96913 20.96913 18.91721 18.26594 18.26594 16.30436 16.30436 18.79606 18.79606 18.79606 19.80604 19.80604 27.06275 27.06275 19.68177 19.68177 18.77782 18.77782 16.99969 16.99969 18.77935 18.77935 18.77935 19.87208 19.87208 26.61473
[63] 26.61473 24.92790 23.95655 23.95655 22.77274 22.77274 23.58385 23.58385 23.58385 24.60241 24.60241 29.78621 29.78621
$df
[1] 69
$residual.scale
[1] 78.62265
diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit
plotSpecs2 <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
#geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
plotObs_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted ASVs")
plotObs_pred
predict(chao1Mod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
480.6067 480.6067 353.4462 353.4462 448.7195 264.4622 264.4622 406.6256 382.2290 680.1825 680.1825 553.0219 553.0219 648.2953 464.0380 464.0380 606.2014 606.2014 756.4058 756.4058 629.2453 629.2453 724.5186 540.2614 540.2614 682.4247 658.0281 658.0281 760.2027 760.2027 633.0422 633.0422
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
728.3155 728.3155 544.0583 544.0583 686.2216 686.2216 721.4327 594.2722 594.2722 689.5455 689.5455 505.2883 505.2883 505.2883 647.4517 647.4517 623.0550 623.0550 624.4071 624.4071 497.2465 497.2465 592.5199 592.5199 408.2626 408.2626 408.2626 550.4260 550.4260 526.0294 526.0294 502.0359
65 66 67 68 69 70 71 72 73 74 75
374.8753 374.8753 470.1487 470.1487 285.8914 285.8914 285.8914 428.0548 428.0548 403.6582 403.6582
$se.fit
[1] 96.26104 96.26104 93.63656 93.63656 93.13686 99.08674 99.08674 104.65049 121.40314 67.52582 67.52582 65.32740 65.32740 61.35005 71.47751 71.47751 76.69512 76.69512 68.67338 68.67338 67.02618 67.02618 61.47610 71.64408 71.64408 75.63246 100.49907 100.49907
[29] 68.66562 68.66562 66.90815 66.90815 60.49784 60.49784 69.94644 69.94644 73.38642 73.38642 66.20526 63.92598 63.92598 57.06097 57.06097 65.78128 65.78128 65.78128 69.31593 69.31593 94.71251 94.71251 68.88102 68.88102 65.71742 65.71742 59.49445 59.49445
[57] 65.72280 65.72280 65.72280 69.54705 69.54705 93.14456 93.14456 87.24108 83.84161 83.84161 79.69860 79.69860 82.53725 82.53725 82.53725 86.10195 86.10195 104.24388 104.24388
$df
[1] 69
$residual.scale
[1] 275.1586
diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit
plotChao1_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred
predict(shanMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
4.524800 4.524800 4.344339 4.344339 4.334988 4.022197 4.022197 4.154170 4.370926 4.974911 4.974911 4.794451 4.794451 4.785100 4.472309 4.472309 4.604282 4.604282 5.165240 5.165240 4.984780 4.984780 4.975429 4.662638 4.662638 4.794611 5.011367 5.011367 5.207335 5.207335 5.026874 5.026874
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
5.017523 5.017523 4.704732 4.704732 4.836705 4.836705 5.152692 4.972232 4.972232 4.962881 4.962881 4.650090 4.650090 4.650090 4.782063 4.782063 4.998819 4.998819 4.985785 4.985785 4.805324 4.805324 4.795974 4.795974 4.483183 4.483183 4.483183 4.615156 4.615156 4.831912 4.831912 4.762373
65 66 67 68 69 70 71 72 73 74 75
4.581912 4.581912 4.572561 4.572561 4.259770 4.259770 4.259770 4.391743 4.391743 4.608499 4.608499
$se.fit
[1] 0.1604232 0.1604232 0.1560494 0.1560494 0.1552166 0.1651324 0.1651324 0.1744046 0.2023236 0.1125347 0.1125347 0.1088710 0.1088710 0.1022425 0.1191204 0.1191204 0.1278158 0.1278158 0.1144472 0.1144472 0.1117021 0.1117021 0.1024526 0.1193980 0.1193980 0.1260448 0.1674861 0.1674861
[29] 0.1144342 0.1144342 0.1115054 0.1115054 0.1008223 0.1008223 0.1165688 0.1165688 0.1223017 0.1223017 0.1103340 0.1065354 0.1065354 0.0950946 0.0950946 0.1096274 0.1096274 0.1096274 0.1155180 0.1155180 0.1578425 0.1578425 0.1147932 0.1147932 0.1095210 0.1095210 0.0991501 0.0991501
[57] 0.1095299 0.1095299 0.1095299 0.1159032 0.1159032 0.1552295 0.1552295 0.1453911 0.1397257 0.1397257 0.1328212 0.1328212 0.1375519 0.1375519 0.1375519 0.1434927 0.1434927 0.1737270 0.1737270
$df
[1] 69
$residual.scale
[1] 0.4585638
diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit
plotShannonH_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred
predict(pielouMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
0.019119982 0.019119982 0.021726336 0.021726336 0.020969819 0.024752823 0.024752823 0.022023135 0.018632849 0.010604187 0.010604187 0.013210541 0.013210541 0.012454023 0.016237028 0.016237028 0.013507339 0.013507339 0.006862322 0.006862322 0.009468676 0.009468676 0.008712158 0.012495163
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
0.012495163 0.009765474 0.006375189 0.006375189 0.005809293 0.005809293 0.008415647 0.008415647 0.007659129 0.007659129 0.011442134 0.011442134 0.008712445 0.008712445 0.006592059 0.009198414 0.009198414 0.008441896 0.008441896 0.012224901 0.012224901 0.012224901 0.009495212 0.009495212
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
0.006104927 0.006104927 0.009352450 0.009352450 0.011958804 0.011958804 0.011202287 0.011202287 0.014985291 0.014985291 0.014985291 0.012255603 0.012255603 0.008865318 0.008865318 0.013176838 0.015783192 0.015783192 0.015026675 0.015026675 0.018809679 0.018809679 0.018809679 0.016079991
73 74 75
0.016079991 0.012689706 0.012689706
$se.fit
[1] 0.005295216 0.005295216 0.005150846 0.005150846 0.005123358 0.005450654 0.005450654 0.005756710 0.006678255 0.003714522 0.003714522 0.003593590 0.003593590 0.003374800 0.003931900 0.003931900 0.004218915 0.004218915 0.003777648 0.003777648 0.003687038 0.003687038 0.003381734
[24] 0.003941063 0.003941063 0.004160460 0.005528345 0.005528345 0.003777221 0.003777221 0.003680545 0.003680545 0.003327921 0.003327921 0.003847678 0.003847678 0.004036908 0.004036908 0.003641880 0.003516499 0.003516499 0.003138862 0.003138862 0.003618557 0.003618557 0.003618557
[47] 0.003812994 0.003812994 0.005210033 0.005210033 0.003789070 0.003789070 0.003615044 0.003615044 0.003272725 0.003272725 0.003615340 0.003615340 0.003615340 0.003825708 0.003825708 0.005123781 0.005123781 0.004799037 0.004612036 0.004612036 0.004384134 0.004384134 0.004540285
[70] 0.004540285 0.004540285 0.004736375 0.004736375 0.005734343 0.005734343
$df
[1] 69
$residual.scale
[1] 0.01513617
diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit
plot_pielouJ_pred <- diversityData %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred
plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).
plotPredictions
ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)
plot_grid(plotObs, plotChao1, plotShan, plotPielou,
plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
nrow = 2, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (geom_segment).
alphaSummary <- tibble(
metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(obsMod, chao1Mod, shanMod, pielouMod)
)
alphaSummary <- alphaSummary %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
alphaSummary <- alphaSummary %>%
select(-model) %>%
unnest(df)
alphaSummary <- alphaSummary %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
bold(i = ~ p< 0.05, j = "p") %>%
colformat_md() %>%
set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
Metric | Term | Estimate | Standard | T Value | p |
Observed ASVs | Intercept | 2.9 x 105 | 1.4 x 105 | 2.12 | 0.038 |
log(Size Class) | 3.2 x 101 | 8.4 x 100 | 3.84 | < 0.001 | |
log(Size Class)2 | -6.4 x 100 | 1.6 x 100 | -4.09 | < 0.001 | |
Latitude | -1.5 x 104 | 7.2 x 103 | -2.12 | 0.038 | |
Latitude2 | 2.0 x 102 | 9.4 x 101 | 2.12 | 0.038 | |
Depth | 5.2 x 100 | 3.3 x 100 | 1.56 | 0.124 | |
Richness (Chao1) | Intercept | 1.1 x 106 | 4.8 x 105 | 2.22 | 0.030 |
log(Size Class) | 8.6 x 101 | 2.9 x 101 | 2.93 | 0.005 | |
log(Size Class)2 | -1.8 x 101 | 5.4 x 100 | -3.31 | 0.001 | |
Latitude | -5.6 x 104 | 2.5 x 104 | -2.22 | 0.030 | |
Latitude2 | 7.3 x 102 | 3.3 x 102 | 2.22 | 0.030 | |
Depth | 2.2 x 101 | 1.2 x 101 | 1.92 | 0.059 | |
Diversity (Shannon H) | Intercept | 1.7 x 103 | 8.1 x 102 | 2.09 | 0.040 |
log(Size Class) | 2.0 x 10-1 | 4.9 x 10-2 | 4.08 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 9.1 x 10-3 | -4.03 | < 0.001 | |
Latitude | -8.7 x 101 | 4.2 x 101 | -2.08 | 0.041 | |
Latitude2 | 1.1 x 100 | 5.5 x 10-1 | 2.08 | 0.041 | |
Depth | 2.1 x 10-2 | 1.9 x 10-2 | 1.09 | 0.281 | |
Evenness (Pielou J) | Intercept | -2.5 x 101 | 2.7 x 101 | -0.94 | 0.349 |
log(Size Class) | -3.8 x 10-3 | 1.6 x 10-3 | -2.37 | 0.021 | |
log(Size Class)2 | 6.6 x 10-4 | 3.0 x 10-4 | 2.21 | 0.030 | |
Latitude | 1.3 x 100 | 1.4 x 100 | 0.94 | 0.349 | |
Latitude2 | -1.7 x 10-2 | 1.8 x 10-2 | -0.94 | 0.350 | |
Depth | -3.7 x 10-4 | 6.4 x 10-4 | -0.57 | 0.569 |
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))
Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.
diversityDataWB <- full_join(diversityData,
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
mutate(pielouJ2 = shannonH/break_estimate) %>%
identity()
diversityDataWB
pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)
Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-0.013926 -0.005094 -0.002507 0.000907 0.105946
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.955e+01 2.753e+01 -0.710 0.4800
log(Size_Class) -3.291e-03 1.668e-03 -1.974 0.0524 .
I(log(Size_Class)^2) 5.747e-04 3.103e-04 1.852 0.0683 .
lat 1.017e+00 1.434e+00 0.709 0.4809
I(lat^2) -1.321e-02 1.867e-02 -0.707 0.4818
depth -2.353e-04 6.605e-04 -0.356 0.7228
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01568 on 69 degrees of freedom
Multiple R-squared: 0.06733, Adjusted R-squared: -0.0002567
F-statistic: 0.9962 on 5 and 69 DF, p-value: 0.4267
Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.
plotBreak <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Richness (Breakaway)")
plotBreak
plotPielou2 <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Evenness (PielouJ)")
plotPielou2
predict(pielouMod2, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0.0116882006 0.0116882006 0.0135742540 0.0135742540 0.0133932226 0.0160282544 0.0160282544 0.0140111437 0.0099141263 0.0043215048 0.0043215048 0.0062075582 0.0062075582 0.0060265267 0.0086615586 0.0086615586 0.0066444478 0.0066444478 0.0010939134 0.0010939134
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
0.0029799668 0.0029799668 0.0027989353 0.0054339672 0.0054339672 0.0034168564 -0.0006801609 -0.0006801609 0.0002000096 0.0002000096 0.0020860630 0.0020860630 0.0019050316 0.0019050316 0.0045400634 0.0045400634 0.0025229527 0.0025229527 0.0008938114 0.0027798648
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
0.0027798648 0.0025988333 0.0025988333 0.0052338652 0.0052338652 0.0052338652 0.0032167544 0.0032167544 -0.0008802629 -0.0008802629 0.0033080946 0.0033080946 0.0051941480 0.0051941480 0.0050131165 0.0050131165 0.0076481484 0.0076481484 0.0076481484 0.0056310377
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
0.0056310377 0.0015340203 0.0015340203 0.0066431363 0.0085291897 0.0085291897 0.0083481582 0.0083481582 0.0109831901 0.0109831901 0.0109831901 0.0089660794 0.0089660794 0.0048690620 0.0048690620
$se.fit
[1] 0.005484918 0.005484918 0.005335376 0.005335376 0.005306903 0.005645925 0.005645925 0.005962946 0.006917505 0.003847596 0.003847596 0.003722331 0.003722331 0.003495703 0.004072762 0.004072762 0.004370059 0.004370059 0.003912983 0.003912983 0.003819127 0.003819127 0.003502885
[24] 0.004082253 0.004082253 0.004309509 0.005726399 0.005726399 0.003912541 0.003912541 0.003812401 0.003812401 0.003447144 0.003447144 0.003985522 0.003985522 0.004181531 0.004181531 0.003772351 0.003642478 0.003642478 0.003251313 0.003251313 0.003748193 0.003748193 0.003748193
[47] 0.003949595 0.003949595 0.005396683 0.005396683 0.003924814 0.003924814 0.003744554 0.003744554 0.003389971 0.003389971 0.003744860 0.003744860 0.003744860 0.003962765 0.003962765 0.005307342 0.005307342 0.004970964 0.004777263 0.004777263 0.004541196 0.004541196 0.004702942
[70] 0.004702942 0.004702942 0.004906057 0.004906057 0.005939777 0.005939777
$df
[1] 69
$residual.scale
[1] 0.01567843
diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit
plot_pielouJ2_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred
plotBreakaway <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) +
scale_y_log10() +
ylab("Richness (Breakaway)")
plotBreakaway
#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame
Why are these not smooth curves?!! What if I redo the model, this time with the same data frame
breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat + I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()
Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-2974.5 -1191.2 -151.6 599.9 6768.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7124615.61 3339862.88 2.133 0.0365 *
log(Size_Class) 244.45 202.35 1.208 0.2312
I(log(Size_Class)^2) -75.16 37.65 -1.996 0.0498 *
lat -370568.38 174035.93 -2.129 0.0368 *
I(lat^2) 4817.28 2265.61 2.126 0.0371 *
depth 151.10 80.15 1.885 0.0636 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared: 0.1414, Adjusted R-squared: 0.0792
F-statistic: 2.273 on 5 and 69 DF, p-value: 0.0567
Note the non statistical significance overall
#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit
plot_break_pred <- diversityDataWB %>%
filter(Depth %in% c("Surface", "Bottom")) %>%
# filter(Station == 4.3) %>%
ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred
plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB
ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)
Summary table I want both breakaway metrics here
bettaTable <- myBet$table %>%
as.data.frame() %>%
rename(estimate = Estimates,
`std.error` = `Standard Errors`,
`p.value`=`p-values`
) %>%
mutate(`statistic` = NA) %>%
rownames_to_column(var = "term") %>%
select(term, estimate, std.error, statistic, p.value) %>%
as_tibble()
bettaTable
alphaSummary2 <- tibble(
metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(breakLm, shanMod, pielouMod2)
)
alphaSummary2 <- alphaSummary2 %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
## Add in willis variables
breakawaySummary <- tibble(
metric = "Richness (Breakaway -- Betta)",
model = NULL,
df = list(bettaTable)
)
alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)
alphaSummary2 <- alphaSummary2 %>%
select(-model) %>%
unnest(df)
alphaSummary2 <- alphaSummary2 %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary2
alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
alphaTable2
Metric | Term | Estimate | Standard | T Value | p |
Richness (Breakaway – Betta) | Intercept | 7.1 x 106 | 2.4 x 102 | NA | < 0.001 |
log(Size Class) | 1.2 x 102 | 6.1 x 101 | NA | 0.058 | |
log(Size Class)2 | -5.0 x 101 | 1.2 x 101 | NA | < 0.001 | |
Latitude | -3.7 x 105 | 6.1 x 100 | NA | < 0.001 | |
Latitude2 | 4.8 x 103 | 1.6 x 10-1 | NA | < 0.001 | |
Depth | 1.5 x 102 | 1.0 x 101 | NA | < 0.001 | |
Richness (Breakaway – LM) | Intercept | 7.1 x 106 | 3.3 x 106 | 2.13 | 0.036 |
log(Size Class) | 2.4 x 102 | 2.0 x 102 | 1.21 | 0.231 | |
log(Size Class)2 | -7.5 x 101 | 3.8 x 101 | -2.00 | 0.050 | |
Latitude | -3.7 x 105 | 1.7 x 105 | -2.13 | 0.037 | |
Latitude2 | 4.8 x 103 | 2.3 x 103 | 2.13 | 0.037 | |
Depth | 1.5 x 102 | 8.0 x 101 | 1.89 | 0.064 | |
Diversity (Shannon H) | Intercept | 1.7 x 103 | 8.1 x 102 | 2.09 | 0.040 |
log(Size Class) | 2.0 x 10-1 | 4.9 x 10-2 | 4.08 | < 0.001 | |
log(Size Class)2 | -3.7 x 10-2 | 9.1 x 10-3 | -4.03 | < 0.001 | |
Latitude | -8.7 x 101 | 4.2 x 101 | -2.08 | 0.041 | |
Latitude2 | 1.1 x 100 | 5.5 x 10-1 | 2.08 | 0.041 | |
Depth | 2.1 x 10-2 | 1.9 x 10-2 | 1.09 | 0.281 | |
Evenness (Pielou J) | Intercept | -2.0 x 101 | 2.8 x 101 | -0.71 | 0.480 |
log(Size Class) | -3.3 x 10-3 | 1.7 x 10-3 | -1.97 | 0.052 | |
log(Size Class)2 | 5.7 x 10-4 | 3.1 x 10-4 | 1.85 | 0.068 | |
Latitude | 1.0 x 100 | 1.4 x 100 | 0.71 | 0.481 | |
Latitude2 | -1.3 x 10-2 | 1.9 x 10-2 | -0.71 | 0.482 | |
Depth | -2.4 x 10-4 | 6.6 x 10-4 | -0.36 | 0.723 |
alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))
myBet$table
plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)